In [1]:
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (10, 6)
In [2]:
from hrv.io import read_from_text
rri = read_from_text("data/08012805.txt")
In [3]:
rri.info()
The RRi series has 2380 values and approximately 20 minutes.
In [4]:
fig, ax = rri.plot()
This is an RRi series recorded during a maximal effort exercise. The first 180 seconds (3 minutes) the subject is resting, after this period the exercise test started and the workload was incremented each minute until the subject's fatigue. Following the exercise, there is a recovery period of approximately 600s (10 minutes).
Looks like there are some noise in the RRi signal. Let's try to filter the time series:
In [5]:
from hrv.filters import quotient, moving_median
In [6]:
fq_rri = quotient(rri)
fig, ax = fq_rri.plot()
In [7]:
fmm_rri = moving_median(rri, order=5)
fig, ax = fmm_rri.plot()
Both filters removed the spikes, but seems that the quotient filter preserved the signal, only removing the noise, while the moving_median filtered the whole tachogram. Let's keep the quotient filter results.
In [8]:
rest_rri = fq_rri.time_range(start=0, end=180)
fig, ax = rest_rri.plot()
In [9]:
from hrv.classical import frequency_domain, time_domain
rest_time_domain = time_domain(rest_rri)
rest_time_domain
Out[9]:
Before extracting Frequency Domain features lets first remove the slow trend from the RRi signal:
In [10]:
from hrv.detrend import polynomial_detrend
detrended_rest_rri = polynomial_detrend(rest_rri, degree=3)
fig, ax = detrended_rest_rri.plot()
Note how the Y-axis is now centered on zero.
In [11]:
detrended_rest_rri.info()
Once our rest signal has only 167 points, lets reduce the segment size and the overlap size of Welch's method to 64 and 32, respectively.
In [12]:
rest_freq_domain = frequency_domain(
detrended_rest_rri,
method="welch",
nperseg=64,
noverlap=32,
interp_method="cubic",
window="hanning",
fs=4.0
)
rest_freq_domain
Out[12]:
In [13]:
recovery_rri = rri.time_range(start=rri.time[-1] - 180, end=rri.time[-1]).reset_time()
fig, ax = recovery_rri.plot()
In [14]:
recovery_rri.info()
In [15]:
recovery_time_domain = time_domain(recovery_rri)
recovery_time_domain
Out[15]:
In [16]:
detrended_recovery_rri = polynomial_detrend(recovery_rri, degree=3)
In [17]:
fig, ax = detrended_recovery_rri.plot()
In [18]:
recovery_freq_domain = frequency_domain(
detrended_recovery_rri,
method="welch",
nperseg=64,
noverlap=32,
interp_method="cubic",
window="hanning",
fs=4.0
)
recovery_freq_domain
Out[18]:
In [19]:
def compare_indices(ax, cond_1, cond_2, index_name, title, y_label):
ax.bar([0, 1], [cond_1[index_name], cond_2[index_name]], color=["b", "r"])
ax.set_xticks([0, 1])
ax.set_xticklabels(["Rest", "Recovery"])
ax.set(ylabel=y_label)
ax.set(title=title)
In [20]:
fig, ax = plt.subplots(2, 2)
fig.set_size_inches(15, 12)
compare_indices(
ax[0][0],
rest_time_domain,
recovery_time_domain,
"rmssd",
title="Time Domain",
y_label="RMSSD (ms)"
)
compare_indices(
ax[0][1],
rest_time_domain,
recovery_time_domain,
"pnn50",
title="Time Domain",
y_label="pNN50 (%)"
)
compare_indices(
ax[1][0],
rest_freq_domain,
recovery_freq_domain,
"hf",
title="Frequency Domain",
y_label="HF (ms²)"
)
compare_indices(
ax[1][1],
rest_freq_domain,
recovery_freq_domain,
"lf",
title="Frequency Domain",
y_label="LF (ms²)"
)
The figure above depicts the comparison between RMSSD, pNN50, HF, and LF extracted on the Rest (blue) and Recovery (red) periods. The reduced values of these indices in the recovery period might indicate that the vagal activity is, at least, partially suppressed after the maximal effort exercise.
The reduced LF (ms²) measure indicates that the RRi series at the recovery period has fewer overall fluctuations compared to the Rest period.
Methods of assessment of the post-exercise cardiac autonomic recovery: A methodological review
Absence of parasympathetic reactivation after maximal exercise
One of the reasons for selecting the Rest and the Recovery periods is due to its stationary behavior. Classical HRV indices expect that the statistical properties of the RRi signal are stable as a function of time. Therefore, extracting classical indices (Time and Frequency domain) in non-stationary segments might bring misleading results.
Let's take a look at the RRi series at the peak of the maximal effort exercise:
In [21]:
peak_exercise_rri = rri.time_range(start=400, end=600)
fig, ax = peak_exercise_rri.plot()
As shown in the above picture, the RRi series during exercise is non-stationary and for this reason, classical analyses are not recommended.
To overcome the non-stationary behavior and also extract information about the dynamics of the HRV in experiments involving physical exercise, Tilt maneuver it is possible to use time-varying method, which consists of splitting the RRi signal into smaller segments (ex: 30s) and calculate the time domain indices of each adjacent segment.
There are also Frequency domain analyses in adjacent smaller segments of the RRi signal like Short Time Fourier Transform, but it is still a work in progress in the hrv module.
In [22]:
from hrv.nonstationary import time_varying
In [23]:
tv_results = time_varying(fq_rri, seg_size=30, overlap=0)
fig, ax = tv_results.plot(index="rmssd", marker="o", color="k")